Voineagu dataset GSE28521

Load data

Cannot use getGEO because everything is in different places, use lumiBatch object to gather all the information

if(!file.exists('./../Data/Voineagu/LumiBatch_Voineagu.RData')){
  
  
  create_Voineagus_pData_df = function(){
    # Output: Dataframe with pData of Voineagu's raw data
   
    # BrainBank, chip_array
    Voineagu_suppMat = read.csv('../Data/Voineagu/aux/supplementary_table_2.csv')
    Voineagu_suppMat$chip_array = paste(Voineagu_suppMat$Chip, Voineagu_suppMat$Array, sep='_')
    Voineagu_suppMat = Voineagu_suppMat[,c('Brain.Bank.Case.Id','Cortex.area','Disease.status',
                                           'chip_array')]
    Voineagu_suppMat$Cortex.area = as.character(Voineagu_suppMat$Cortex.area)
    
    # MATCH ALL FRONTAL + TEMPORAL IDS:
    # GSM, BrainBank (relation obtained from the series matrix phenotype information)
    ids_mapping = read.delim('../Data/Voineagu/aux/ids_mapping.txt', header=FALSE)
    colnames(ids_mapping)[1:2] = c('GSM','Brain.Bank.Plus')
    id_split = strsplit(as.character(ids_mapping$Brain.Bank.Plus),'_')
    ids_mapping$Disease.status = sapply(id_split, `[`, 1)
    ids_mapping$Brain.Bank.Case.Id = sapply(id_split, `[`, 2)
    ids_mapping$Cortex.area = sapply(id_split, `[`, 3)
    
    cortex_area = list('C'='cerebellum','T'='temporal','F'='frontal')
    ids_mapping$Cortex.area[ids_mapping$Cortex.area %in% names(cortex_area)] = 
      cortex_area[ids_mapping$Cortex.area[ids_mapping$Cortex.area %in% names(cortex_area)]]
    
    disease_status = list('A'='autism', 'C'='control')
    ids_mapping$Disease.status[ids_mapping$Disease.status %in% names(disease_status)] = 
      disease_status[ids_mapping$Disease.status[ids_mapping$Disease.status %in% names(disease_status)]]
    
    # BrainBank, metadata
    samples_metadata = read.delim2('../Data/Voineagu/aux/supplementary_table_1.tsv')
    
    # Frontal and temporal info:
    ids_mapping_f_t = merge(Voineagu_suppMat, ids_mapping, 
                            by=c('Brain.Bank.Case.Id','Disease.status','Cortex.area'))
    
    # MATCH ALL CEREBELLUM IDS:
    # GSM, BrainBank, chip_array (cerebellum)
    ids_mapping_c = read.delim('../Data/Voineagu/aux/cerebellum_ids_mapping.txt', header=FALSE)
    colnames(ids_mapping_c) = c('GSM','Brain.Bank.Plus','chip_array')
    
    id_split = strsplit(as.character(ids_mapping_c$Brain.Bank.Plus),'_')
    ids_mapping_c$Disease.status = sapply(id_split, `[`, 1)
    ids_mapping_c$Brain.Bank.Case.Id = sapply(id_split, `[`, 2)
    ids_mapping_c$Cortex.area = sapply(id_split, `[`, 3)
    
    ids_mapping_c$Cortex.area[ids_mapping_c$Cortex.area %in% names(cortex_area)] = 
      cortex_area[ids_mapping_c$Cortex.area[ids_mapping_c$Cortex.area %in% names(cortex_area)]]
    
    ids_mapping_c$Disease.status[ids_mapping_c$Disease.status %in% names(disease_status)] = 
    disease_status[ids_mapping_c$Disease.status[ids_mapping_c$Disease.status %in% names(disease_status)]]
    
    # JOIN FRONTAL+TEMPORAL INFO WITH CEREBELLUM AND ADD METADATA:
    ids_mapping_all = rbind(ids_mapping_f_t, ids_mapping_c)
    ids_mapping_all = ids_mapping_all[,c('GSM','chip_array','Brain.Bank.Case.Id','Cortex.area',
                                         'Disease.status')]
    
    samples_full_data = merge(ids_mapping_all, samples_metadata, by='Brain.Bank.Case.Id')
    samples_full_data$Cortex.area = as.character(samples_full_data$Cortex.area)
    
    # Add age group data
    age_group = cut(as.numeric(samples_full_data$AGE), c(-0.5, 0, 0.6, 10, 20, 30, 40, 50, 60, 70, 80),
                    labels=c('Fetal','Infant','Child','10-20','20s','30s','40s','50s','60s','70s'))
    samples_full_data$age_group = as.character(age_group)
    
    return(samples_full_data)
  }
  
  # AssayData data
  # Supplementary file GSE28521_non-normalized_data.txt.gz
  LumiBatch_v = lumiR('../Data/Voineagu/GSE28521_non-normalized_data.txt.gz')
  
  # Phenotype data
  phenoData = create_Voineagus_pData_df()
  phenoData = phenoData[match(colnames(exprs(LumiBatch_v)), phenoData$chip_array),]
  phenoData$SEX[phenoData$SEX=='5'] = 'F'
  rownames(phenoData) = phenoData$GSM
  
  # Change subject IDs from chip_array to GSM format
  chipArray_to_gsm = phenoData$GSM
  names(chipArray_to_gsm) = phenoData$chip_array
  #colnames(LumiBatch_v) = chipArray_to_gsm[colnames(LumiBatch_v)]
  
  # LumiBatch object
  pData(LumiBatch_v) = phenoData
  
  save(LumiBatch_v, file='../Data/Voineagu/LumiBatch_Voineagu.RData')
  
  
  remove(assayData, phenoData, chipArray_to_gsm)
} else {
  load(paste('../Data/Voineagu/LumiBatch_Voineagu.RData', sep='/'))
  datExpr = exprs(LumiBatch_v)
  datMeta = pData(LumiBatch_v)
}

# Get Biomart information
ensembl = useMart("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl", host="feb2014.archive.ensembl.org") 
getinfo = c("illumina_humanref_8_v3", "ensembl_gene_id","external_gene_id", "entrezgene", 
            "chromosome_name", "start_position", "end_position")
geneDat = getBM(attributes = getinfo, filters="illumina_humanref_8_v3",
                values = rownames(datExpr), mart = ensembl)
idx = match(rownames(datExpr), geneDat[,"illumina_humanref_8_v3"])
datGenes = cbind(rownames(datExpr), geneDat[idx,])
rownames(datGenes) = datGenes[,1]
datGenes$length = datGenes$end_position-datGenes$start_position

# SFARI Genes
SFARI_genes = read_csv('./../Data/SFARI/SFARI_genes_with_ensembl_IDs.csv')


# GO Annotations
GO_annotations = read.csv('./../Data/GO_annotations/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID' = as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal' = 1)


# Combine SFARI and GO information
gene_info = data.frame('ID'=datGenes$ensembl_gene_id) %>% left_join(SFARI_genes, by='ID') %>% 
            mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
            left_join(GO_neuronal, by='ID') %>% 
            mutate('gene.score'=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`)) %>%
            mutate('gene.score'=ifelse(is.na(gene.score), 'None', gene.score))

rm(ensembl, geneDat, GO_annotations, LumiBatch_v, getinfo, idx, create_Voineagus_pData_df)


Remove cerebellum samples:

datMeta = datMeta[datMeta$Cortex.area!='cerebellum',]
datExpr = datExpr[,datMeta$chip_array]

Check sample distribution

RNA-Seq for 58 cortical brain-tissue samples across frontal and temporal lobes, comprising 29 samples from control subjects and 29 samples from ASD subjects

print(paste0('Dataset includes ', nrow(datExpr), ' genes from ', ncol(datExpr), ' samples belonging to ', length(unique(datMeta$Brain.Bank.Case.Id)), ' different subjects.'))
## [1] "Dataset includes 24526 genes from 58 samples belonging to 32 different subjects."


Sex distribution: There are five times more males than females

table(datMeta$SEX)
## 
##  F  M 
##  9 49


Age distribution: Subjects between 5 and 56 years old with a mean close to 29

summary(datMeta$AGE)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00   18.25   28.50   28.79   39.00   56.00


Diagnosis distribution: Fairly balanced

table(datMeta$Disease.status)
## 
##  autism control 
##      29      29


Brain region distribution:

table(datMeta$Cortex.area)
## 
##  frontal temporal 
##       32       26


Diagnosis and brain region are very well balanced

table(datMeta$Disease.status, datMeta$Cortex.area)
##          
##           frontal temporal
##   autism       16       13
##   control      16       13


Filtering


1. Filter genes with start or end position missing

to_keep = !is.na(datGenes$length)
datGenes = datGenes[to_keep,]
datExpr = datExpr[to_keep,]
#rownames(datGenes) = datGenes$ensembl_gene_id

print(paste0('Removed ', sum(!to_keep), ' genes, ', sum(to_keep), ' remaining'))
## [1] "Removed 780 genes, 23746 remaining"


2. Filter genes with low expression levels

The density distribution doesn’t seem to suggest any clearly defined threshold. Selecting 300

plot_data = data.frame('id'=rownames(datExpr), 'mean_expression' = rowMeans(datExpr))
ggplotly(plot_data %>% ggplot(aes(x=mean_expression)) + 
         geom_density(color='#0099cc', fill='#0099cc', alpha=0.3) + 
         geom_vline(xintercept=300, color='gray') + scale_x_log10() + 
         ggtitle('gene Mean Expression distribution') + theme_minimal())
to_keep = rowMeans(datExpr)>300
datGenes = datGenes[to_keep,]
datExpr = datExpr[to_keep,]

print(paste0('Removed ', sum(!to_keep), ' genes, ', sum(to_keep), ' remaining'))
## [1] "Removed 8943 genes, 14803 remaining"


  1. Filter outlier samples

Using node connectivity as a distance measure, normalising it and filtering out genes farther away than 2 standard deviations from the left (lower connectivity than average, not higher)

absadj = datExpr %>% bicor %>% abs
## alpha: 1.000000
netsummary = fundamentalNetworkConcepts(absadj)
ku = netsummary$Connectivity
z.ku = (ku-mean(ku))/sqrt(var(ku))

plot_data = data.frame('ID'=rownames(datMeta), 'sample'=1:length(z.ku), 'distance'=z.ku, 'Region'=datMeta$Cortex.area, 'Diagnosis'=datMeta$Disease.status)
ggplotly(plot_data %>% ggplot(aes(sample, distance, color=Diagnosis)) + geom_point() + geom_hline(yintercept=-2, color='gray') + theme_minimal())
print(paste0('Outlier samples: ', paste(as.character(plot_data$ID[plot_data$distance< -2]), collapse=', ')))
## [1] "Outlier samples: GSM706448, GSM706425, GSM706447, GSM706455"
to_keep = abs(z.ku)<2
datMeta = datMeta[to_keep,]
datExpr = datExpr[,to_keep]

print(paste0('Removed ', sum(!to_keep), ' samples, ', sum(to_keep), ' remaining'))
## [1] "Removed 4 samples, 54 remaining"
rm(absadj, netsummary, ku, z.ku, plot_data, to_keep)
print(paste0('After filtering, the dataset consists of ', nrow(datExpr), ' genes and ', ncol(datExpr), ' samples'))
## [1] "After filtering, the dataset consists of 14803 genes and 54 samples"



Normalisation

# Transform Signal to Ratio
#row_mean = apply(datExpr, 1, mean)
#ratio = (datExpr)/(row_mean+1)
#exprs(LumiBatch_v) = as.matrix(ratio[,colnames(ratio) %in% colnames(datExpr)])
#exprs = as.matrix(ratio[,colnames(ratio) %in% colnames(datExpr)])
#LumiBatch_v = ExpressionSet(assayData = exprs)

LumiBatch_v = ExpressionSet(assayData = datExpr)
rownames(datMeta) = datMeta$chip_array
pData(LumiBatch_v) = datMeta

# Variance Stabilisation
LumiBatch_v = lumiT(LumiBatch_v, method = 'log2', ifPlot = TRUE)
## Perform log2 transformation ...
# Normalisation
LumiBatch_v = lumiN(LumiBatch_v, method='rsn')
## Perform rsn normalization ...
## 2019-09-03 23:13:45 , processing array  1 
## 2019-09-03 23:13:45 , processing array  2 
## 2019-09-03 23:13:46 , processing array  3 
## 2019-09-03 23:13:46 , processing array  4 
## 2019-09-03 23:13:46 , processing array  5 
## 2019-09-03 23:13:46 , processing array  6 
## 2019-09-03 23:13:46 , processing array  7 
## 2019-09-03 23:13:46 , processing array  8 
## 2019-09-03 23:13:47 , processing array  9 
## 2019-09-03 23:13:47 , processing array  10 
## 2019-09-03 23:13:47 , processing array  11 
## 2019-09-03 23:13:47 , processing array  12 
## 2019-09-03 23:13:47 , processing array  13 
## 2019-09-03 23:13:47 , processing array  14 
## 2019-09-03 23:13:47 , processing array  15 
## 2019-09-03 23:13:47 , processing array  16 
## 2019-09-03 23:13:47 , processing array  17 
## 2019-09-03 23:13:47 , processing array  18 
## 2019-09-03 23:13:47 , processing array  19 
## 2019-09-03 23:13:47 , processing array  20 
## 2019-09-03 23:13:47 , processing array  21 
## 2019-09-03 23:13:47 , processing array  22 
## 2019-09-03 23:13:47 , processing array  23 
## 2019-09-03 23:13:48 , processing array  24 
## 2019-09-03 23:13:48 , processing array  25 
## 2019-09-03 23:13:48 , processing array  26 
## 2019-09-03 23:13:48 , processing array  27 
## 2019-09-03 23:13:48 , processing array  28 
## 2019-09-03 23:13:48 , processing array  29 
## 2019-09-03 23:13:48 , processing array  30 
## 2019-09-03 23:13:48 , processing array  31 
## 2019-09-03 23:13:48 , processing array  32 
## 2019-09-03 23:13:48 , processing array  33 
## 2019-09-03 23:13:48 , processing array  34 
## 2019-09-03 23:13:48 , processing array  35 
## 2019-09-03 23:13:48 , processing array  36 
## 2019-09-03 23:13:48 , processing array  37 
## 2019-09-03 23:13:49 , processing array  38 
## 2019-09-03 23:13:49 , processing array  39 
## 2019-09-03 23:13:49 , processing array  40 
## 2019-09-03 23:13:49 , processing array  41 
## 2019-09-03 23:13:49 , processing array  42 
## 2019-09-03 23:13:49 , processing array  43 
## 2019-09-03 23:13:49 , processing array  44 
## 2019-09-03 23:13:49 , processing array  45 
## 2019-09-03 23:13:49 , processing array  46 
## 2019-09-03 23:13:49 , processing array  47 
## 2019-09-03 23:13:49 , processing array  48 
## 2019-09-03 23:13:49 , processing array  49 
## 2019-09-03 23:13:49 , processing array  50 
## 2019-09-03 23:13:49 , processing array  51 
## 2019-09-03 23:13:49 , processing array  52 
## 2019-09-03 23:13:50 , processing array  53 
## 2019-09-03 23:13:50 , processing array  54
save(LumiBatch_v, file=paste0('../Data/Voineagu/LumiBatch_v_preprocessed.RData'))
# Update datExpr rownames
rownames(datExpr) = datGenes[rownames(datExpr),]$ensembl_gene_id

datExpr = datExpr[!is.na(rownames(datExpr)),]

datExpr = log2(datExpr)

DEA

mod = model.matrix(~Disease.status, data=datMeta)
corfit = duplicateCorrelation(datExpr, mod, block=datMeta$Brain.Bank.Case.Id)
fit = eBayes(lmFit(datExpr, mod,block=datMeta$Brain.Bank.Case.Id, 
                   correlation = corfit$consensus),trend=T)
DE = topTable(fit,coef=2, number=Inf, sort.by="none")

plot_data = DE %>% left_join(gene_info, by='ID')

ggplotly(plot_data[abs(plot_data$logFC)<1,] %>% ggplot(aes(x=gene.score, y=abs(logFC), fill=gene.score)) + 
         geom_boxplot() + scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + 
         theme_minimal() + theme(legend.position='none') + xlab('Gene score') + ylab('abs(lfc)'))

Mean and SD by SFARI score

No clear relation between SFARI score and mean expression/SD. This dataset doesn’t seem to have the same problem as Gandal’s

plot_data = data.frame('ID'=rownames(datExpr), 'MeanExpr'=rowMeans(datExpr), 
                       'SDExpr'=apply(datExpr,1,sd)) %>% 
            left_join(gene_info, by='ID')

p1 = ggplotly(plot_data %>% ggplot(aes(gene.score, MeanExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
              theme(legend.position='none') + xlab('Gene score') + ylab('mean expression'))

p2 = ggplotly(plot_data %>% ggplot(aes(gene.score, SDExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
              ggtitle('Mean Expression (left) and SD (right) by SFARI score') + 
              theme(legend.position='none') + xlab('Gene score') + ylab('standard deviation'))

subplot(p1, p2, nrows=1)
rm(plot_data, p1, p2)



Data seems to be close to homoscedastic

plot_data = data.frame('ID'=rownames(datExpr), 'Mean'=rowMeans(datExpr), 'SD'=apply(datExpr,1,sd))

plot_data %>% ggplot(aes(Mean, SD)) + geom_point(color='#0099cc', alpha=0.1) + 
  geom_smooth(method='lm', color='gray', size=0.5) + theme_minimal()

rm(plot_data)

Batch Correction

No Batch information available in the metadata!

Visualisations


Samples

PCA: There doesn’t seem to be a specific pattern related to diagnosis or brain region

*There weren’t any patterns in MDS and t-SNE either (MDS was the same as PCA)

pca = datExpr %>% t %>% prcomp

datMeta$ID = datMeta$chip_array

plot_data = data.frame('ID'=colnames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>% 
            left_join(datMeta, by='ID') %>% 
            dplyr::select('PC1','PC2','AGE','age_group','SEX','Disease.status', 'Cortex.area')

p1 = plot_data %>% ggplot(aes(PC1, PC2, color=Disease.status)) + geom_point() + theme_minimal() + 
  ggtitle('PCA') +
              xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
              ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)'))

p2 = plot_data %>% ggplot(aes(PC1, PC2, color=Cortex.area)) + geom_point() + theme_minimal() + ggtitle('PCA') +
              xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) + 
              ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)'))

grid.arrange(p1, p2, nrow=1)

rm(pca, plot_data)


Genes

  • First Principal Component explains 96.6% of the total variance

  • There’s a really strong correlation between the mean expression of a gene and the 1st principal component

pca = datExpr %>% prcomp

plot_data = data.frame( 'PC1' = pca$x[,1], 'PC2' = pca$x[,2], 'MeanExpr'=rowMeans(datExpr))

plot_data %>% ggplot(aes(PC1, PC2, color=MeanExpr)) + geom_point(alpha=0.5) + theme_minimal() + 
              scale_color_viridis() + ggtitle('PCA') +
              xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
              ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)'))

rm(pca, plot_data)



Save preprocessed dataset

save(datExpr, datMeta, datGenes, file='./../Data/Voineagu/preprocessed_data.RData')
#load('./../Data/Voineagu/preprocessed_data.RData')



Session info

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] viridis_0.5.1         viridisLite_0.3.0     gridExtra_2.3        
##  [4] plotly_4.9.0          forcats_0.4.0         stringr_1.4.0        
##  [7] dplyr_0.8.2           purrr_0.3.2           readr_1.3.1          
## [10] tidyr_0.8.3           tibble_2.1.3          ggplot2_3.2.0        
## [13] tidyverse_1.2.1       WGCNA_1.68            fastcluster_1.1.25   
## [16] dynamicTreeCut_1.63-1 limma_3.40.6          biomaRt_2.40.4       
## [19] lumi_2.36.0           GEOquery_2.52.0       Biobase_2.44.0       
## [22] BiocGenerics_0.30.0  
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_0.2.5            robust_0.4-18.1            
##   [3] RSQLite_2.1.2               AnnotationDbi_1.46.1       
##   [5] htmlwidgets_1.3             grid_3.6.1                 
##   [7] BiocParallel_1.18.1         munsell_0.5.0              
##   [9] codetools_0.2-16            preprocessCore_1.46.0      
##  [11] statmod_1.4.32              nleqslv_3.3.2              
##  [13] withr_2.1.2                 colorspace_1.4-1           
##  [15] knitr_1.23                  rstudioapi_0.10            
##  [17] stats4_3.6.1                robustbase_0.93-5          
##  [19] labeling_0.3                GenomeInfoDbData_1.2.1     
##  [21] bit64_0.9-7                 rhdf5_2.28.0               
##  [23] vctrs_0.2.0                 generics_0.0.2             
##  [25] xfun_0.8                    R6_2.4.0                   
##  [27] doParallel_1.0.15           GenomeInfoDb_1.20.0        
##  [29] illuminaio_0.26.0           locfit_1.5-9.1             
##  [31] bitops_1.0-6                reshape_0.8.8              
##  [33] DelayedArray_0.10.0         assertthat_0.2.1           
##  [35] promises_1.0.1              scales_1.0.0               
##  [37] nnet_7.3-12                 gtable_0.3.0               
##  [39] affy_1.62.0                 methylumi_2.30.0           
##  [41] rlang_0.4.0                 zeallot_0.1.0              
##  [43] genefilter_1.66.0           splines_3.6.1              
##  [45] rtracklayer_1.44.3          lazyeval_0.2.2             
##  [47] acepack_1.4.1               impute_1.58.0              
##  [49] broom_0.5.2                 checkmate_1.9.4            
##  [51] BiocManager_1.30.4          yaml_2.2.0                 
##  [53] modelr_0.1.5                crosstalk_1.0.0            
##  [55] GenomicFeatures_1.36.4      backports_1.1.4            
##  [57] httpuv_1.5.1                Hmisc_4.2-0                
##  [59] tools_3.6.1                 nor1mix_1.3-0              
##  [61] affyio_1.54.0               RColorBrewer_1.1-2         
##  [63] siggenes_1.58.0             Rcpp_1.0.1                 
##  [65] plyr_1.8.4                  base64enc_0.1-3            
##  [67] progress_1.2.2              zlibbioc_1.30.0            
##  [69] RCurl_1.95-4.12             prettyunits_1.0.2          
##  [71] rpart_4.1-15                openssl_1.4                
##  [73] bumphunter_1.26.0           S4Vectors_0.22.0           
##  [75] SummarizedExperiment_1.14.1 haven_2.1.1                
##  [77] cluster_2.1.0               magrittr_1.5               
##  [79] data.table_1.12.2           mvtnorm_1.0-11             
##  [81] matrixStats_0.54.0          mime_0.7                   
##  [83] hms_0.5.1                   evaluate_0.14              
##  [85] xtable_1.8-4                XML_3.98-1.20              
##  [87] mclust_5.4.5                readxl_1.3.1               
##  [89] IRanges_2.18.2              compiler_3.6.1             
##  [91] minfi_1.30.0                KernSmooth_2.23-15         
##  [93] crayon_1.3.4                htmltools_0.3.6            
##  [95] later_0.8.0                 mgcv_1.8-28                
##  [97] pcaPP_1.9-73                Formula_1.2-3              
##  [99] rrcov_1.4-7                 lubridate_1.7.4            
## [101] DBI_1.0.0                   MASS_7.3-51.4              
## [103] Matrix_1.2-17               cli_1.1.0                  
## [105] quadprog_1.5-7              GenomicRanges_1.36.0       
## [107] pkgconfig_2.0.2             fit.models_0.5-14          
## [109] GenomicAlignments_1.20.1    registry_0.5-1             
## [111] foreign_0.8-72              xml2_1.2.2                 
## [113] foreach_1.4.7               annotate_1.62.0            
## [115] rngtools_1.4                pkgmaker_0.27              
## [117] multtest_2.40.0             beanplot_1.2               
## [119] XVector_0.24.0              bibtex_0.4.2               
## [121] rvest_0.3.4                 doRNG_1.7.1                
## [123] scrime_1.3.5                digest_0.6.19              
## [125] Biostrings_2.52.0           rmarkdown_1.13             
## [127] base64_2.0                  cellranger_1.1.0           
## [129] htmlTable_1.13.1            DelayedMatrixStats_1.6.0   
## [131] curl_3.3                    shiny_1.3.2                
## [133] Rsamtools_2.0.0             nlme_3.1-141               
## [135] jsonlite_1.6                Rhdf5lib_1.6.0             
## [137] askpass_1.1                 pillar_1.4.2               
## [139] lattice_0.20-38             httr_1.4.0                 
## [141] DEoptimR_1.0-8              survival_2.44-1.1          
## [143] GO.db_3.8.2                 glue_1.3.1                 
## [145] iterators_1.0.12            bit_1.1-14                 
## [147] stringi_1.4.3               HDF5Array_1.12.2           
## [149] blob_1.2.0                  latticeExtra_0.6-28        
## [151] memoise_1.1.0